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Energy injection into the early universe can induce turbulent motions of the primordial plasma, 
which in turn act as a source for gravitational radiation. Earlier work computed the amplitude 
and characteristic frequency of the relic gravitational wave background, as a function of the total 

■ energy injected and the stirring scale of the turbulence. This paper computes the frequency spec- 
' trum of relic gravitational radiation from a turbulent source of the stationary Kolmogoroff form 
, which acts for a given duration, making no other approximations. We also show that the limit of 

■ long source wavelengths, commonly employed in aeroacoustic problems, is an excellent approxima- 
5^ ' tion. The gravitational waves from cosmological turbulence around the electroweak energy scale 
CL(| will be detectable by future space-based laser interferometers for a substantial range of turbulence 

. parameters. 

Q\ ■ PACS numbers: 98.70.Vc, 98.80.-k 

(N 



6 



X 
S3 



I. INTRODUCTION 



Direct detection of a relic gravitational wave background is a subject of considerable current interest (see [11,0] i and 
01 Hj 01 for recent reviews) , motivated by planned satellite detection missions in the near future @ . Gravitational 
wave detection could probe directly the physical conditions in the early universe at the epoch of radiation genera- 
tion 0, since after been generated, gravitational radiation freely propagates throughout the entire evolution of the 
\ universe. Once generated, any gravitational wave spectrum retains its shape, with all wavelengths simply scaling 
• with the expansion of the universe. Various possibilities for early-universe physics leading to detectable cosmological 
I gravitational wave backgrounds include quantum fluctuations during inflation Q and subsequent oscillating classi- 
f — , cal fields during reheating cosmological defects bubble wall motions and collisions during phase transitions 
^ ■ [O, [H, d ; plasma turbulence [H, El 111, ^ ; and cosmological magnetic fields [Uli^. Depending on wavelength, 
(y-^ , the resulting gravitational waves might be detected either directly or through their imprint on the polarization of the 
' cosmic microwave background [l^ . If detected, gravitational radiation generated in the earl y u niverse would provide 
t — ] a remarkable new window into physics beyond the standard model of particle physics (e.g., |l9l. [2Qi]). 

■ In this paper we revisit the generation of a cosmological gravitational wave background from turbulent motion of the 
^ ' primordial plasma. We employ methods similar to those originally developed in aeroacoustics for calculating sound 
generation by turbulent flows [21, 22, 23, 24]. This allows us to incorporate the influence of the temporal characteristics 
of turbulent fluctuations on the gravitational wave generation process, and thus to determine the spectrum of the 
emitted gravitational waves at all frequencies. For simplicity, we assume isotropic non-helical turbulence, ignoring 
all possibilities for generating polarized gravitational wave s 12511 . (Polarized radiation might be generated through 
anisotropic stress of the helical primordial rnagnetic field [26[, or from other parity- violating sources in the early 
universe such as Chern-Simons coupling [27l . l28l [29| or an axion field coupling with gravity 30] . Detection of these 
polarized backgrounds are discussed in Ref. |3]|.) 

As is well known, gravitational waves are sourced by the transverse and traceless part of the stress-energy tensor 
(see, e.g., (32j). In our case the stress-energy tensor results from turbulent plasma motions: 

Tij{x) (xwviix)vj{x), (1) 

where v(x) is the velocity vector field of the fluid and w = p + e is the enthalpy density with p and e the pressure 
and the energy density of plasma, which is assumed to be constant throughout space To model a period of 
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cosmological turbulence, we assume that at time in the early universe, a vacuum energy density pvac is converted 
into (turbulent) kinetic energy of the cosmological plasma via stirring on a characteristic source length scale L5, over 
a time scale Tstir [13 ■ The characteristic length scale Ls of the generated fluctuations is directly related to the Hubble 
length = H~^{t^,) at the time of energy injection. We consider only a forward cascade: after being generated on 
the length scale L5, the turbulence kinetic energy cascades from larger to smaller scales. The cascade stops at some 
damping scale Ld, when the energy of turbulence thermalizes due to some dissipation mechanism, such as viscosity 
or plasma resistivity. In this paper we consider w, pvac, '''stir, Ls, and Ld as phenomenological parameters which 
can approximately describe any period of cosmological turbulence, and derive the dependence of the gravitational 
wave spectrum on these parameters. As expected from the universal nature of turbulence, the shape of the spectrum 
scales with the characteristic amplitude and frequency of the gravitational radiation. 

We perform the computation of the gravitational wave spectrum in real space, instead of using conventional Fourier 
space techniques as in Ref. This makes the physical interpretation of all quantities straightforward. The spatial 
structure of the turbulence is taken to be isotropic with a Kolmogoroff spectrum (ssj , and the time dependence of 
the turbulence is described by the Kraichnan time auto-correlation function [s^]- While relativistic turbulence in 
the early universe might depart somewhat from these scalings, these assumptions are based on observed properties of 
laboratory turbulence and will give the correct qualitative features of the resulting radiation spectrum. Generalization 
to alternative turbulence models is straightforward. We use natural units h — c — kB = ^ throughout. 



II. GENERAL FORMALISM 

We assume the duration of the turbulence, tt, is much less than the Hubble time H^^ so the effects of the 

expansion of the universe may be neglected in the generation of gravitational radiation. This adiabatic assumption 
will be valid for any turbulence which is produced in a realistic cosmological phase transition [s^. (Note that the 
duration of the turbulence tt can be substantially longer than the stirring time tsth] see the detailed discussion in 
Then the radiation equation in real space can be written as [s^ 

V2/i,,(x,i) - ^/i«,(x,t) = -167rG5y(x,i). (2) 
where /iy(x, i) is the tensor metric perturbation, the traceless part of the stress-energy tensor Tij{jx.,t) is [s^ 

5.,(x,i) = r,,(x,i)-i<5,,r,!=(x,t), (3) 

and t is physical time. During the period of turbulence, the stress tensor takes the form of Eq. ([T]) (33l |. 
The general solution of Eq. Q is [H ill , 

= agJ d3,'M5:^_|i^. (4) 

Due to the stochastic character of the turbulent stress tensor 5"^ , the generated metric perturbations hij also are 
stochastic. We aim to derive the energy density spectrum of these perturbations at the end of the turbulent phase; 
after that the amplitude and wavelength of the gravitational radiation scales simply with the expansion of the universe. 
The energy density of gravitational waves is defined as [l*] 

PGw{^,t) = _l^(a,/^,,(x,i)9,/.,,(x,t)) = d^xM^x-i^^M^^l^l^^,^^ 

where the brackets denote an ensemble average over realizations of the stochastic source, t' = i — |x — x'| and 
t" = i - |x-x"|. 



A. Localized Source 



We will first consider turbulence in a bounded region of space centered around x = 0. In this case, the energy 
density flux P(x, <) of the radiation propagating outward in the direction n is just 



P(x) = np(x,i)- 



(6) 



At large distances from the turbulent source, the far-field approximation is justified |33.[36|. This assumption replaces 
|x — x'l by |x| in Eq. ([5]), yielding for the gravitational wave energy density flux 

P(x) = j d3x'd3x"(ai5,,(x', t')dtS,,{^\ t")). (7) 

The flux from a spatially bounded source drops as the inverse square of the distance from the radiation source, as 
expected. 

The autocorrelation function of the tensor metric perturbations is defined as 

i(x,T) = ^^{dth,,{^,t)dth,j{^,t + T)), (8) 

with T ~ t' — t, such that pgvf(x) = ^(x, 0). Defining the usual Fourier transform of L(jx., r) as 

/(x,c.) = i-|dre--L(x,r), (9) 
with Lo as the angular frequency, it readily follows that 

PGw{^) ^ J du}I{x,uj), (10) 



(11) 



and therefore /(x, w) represents the spectral energy density of induced gravitational waves [l|,|3 
Substituting Eq. (g]) into Eq. (HJ gives 



l(x,t) = / d3x'dV(a,5,,(x',t')at5.,(x",<")). 

27r|x|-" J 



For the case of stationary turbulence, it can be proven that [23[ 

{^tS^J{K',t')^tS^,{x",t")) - -a2(5y(x',<')^y(x",i")>. (12) 

Using Eq. IT^ with the far-field approximation |x — x'| = |x| — x-x'/|x|, and using the fact that the cross-correlation 
of a stationary random function is independent of time translation, Eq. (jlip reduces to 

i(x,r) = ^^dfj dVd^x"(5.,(x',i)S., (x",r')), (13) 

where 

r' = t + r+i^-(x"-x'). (14) 
|x| 

Defining the two-point time-delayed forth order correlation tensor by 

i?,,feKx',^,T) = ■\{S,,{x',t)SM{^",t + T)), (15) 

w 

where ^ — x" — x' and w = p + p is the enthalpy density of the plasma, Eq. (fT3)) yields 

i(x, r) = ^^9^ J dSx'd^^i?,,,, (^x', I, r + ^ . . (16) 

Fourier transforming this equation gives 

H^, ^) = p j d^x'i/,,,, [^', —uj, ujj (17) 

(summation on i and j assumed), where the four-dimensional power spectral energy density tensor of stationary 
turbulence is defined as 



iJ,,fcKx',k,c.) = ^ j d3|dre'(---'' «)i?,,fei(x',tr). 



(18) 



Equation (fT7|) allows us to calculate the spectral energy density of gravitational waves from a localized source, if the 
real-space statistical properties of the turbulent source are known. 
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B. Spatially Homogeneous Source of Finite Duration 



For a cosmological source of stochastic gravitational radiation, we assume that the source is statistically homoge- 
neous, so that the averaged correlators of the stress tensor have no spatial dependence, and isotropic, so that the 
correlator between two spatial points depends only on the distance between the points and not on the direction. We 
can also simply account for the expansion of the universe by a simple rescaling of the frequency of all radiation after 
its production, so we compute the radiation spectrum in a non-expanding spacetime and include the expansion effect 
at the end. 

With these assumptions, Eq. simplifies to 

i/y7ci(x',k, w) = Hijki{h,u!) 

= / rf'|dTe*--e-* «%fc,($,T) (19) 

= ^/ dTda'e'^^jomR^M^:r), (20) 

so Hijij{'kuj,uj) = Hijij{uj,uj) independent of the observation direction x, as expected on physical grounds. Now 
consider a stochastic source lasting for a finite duration tt, the duration of the turbulent source. The total radiation 
energy spectrum at some point and time is obtained by integrating over all source regions with a light-like separation 
from the observer, which comprises a spherical shell around the observer with a thickness corrresponding to the 
duration of the phase transition, and a radius equal to the proper distance along any light-like path from the observer 
to the source. Due to statistical isotropy and homogeneity, the integral is trivial, contributing only a volume factor, 
giving for the total energy spectrum 

PGw{^) = '^l'^^ = l&ir'^uj^Gw'^TTH^jiAuj.uj). (21) 
a Inw 

This spectrum is of course independent of the position of the observer, as it should be for a stochastic background. 
In the absence of the expansion of the universe, a stochastic source generates a spectrum of radiation which then 
remains constant for all later times. 



III. STATISTICS OF STATIONARY KOLMOGOROFF TURBULENCE 



For a particular model of turbulent motion, the correlations needed for computing gravitational radiation can be 
estimated. Here we consider the simplest turbulence model, the original Kolmogoroff picture. The spectral function 
Fij (k, r) for stationary, isotropic and homogenous turbulence is defined as a spatial Fourier transform of the two-point 
velocity correlation function 

i?,j(r,T) = (v,(x,t)w,(x + r,i + T)). (22) 
This function can be expressed in the form [l^l 



F,,.(k,r) = -4- U,^. - ^ /(,,,,r), (23) 



where is the one-dimensional turbulent spectrum of energy density, rjk is the autocorrelation function 34], and 
the function firik^r) characterizes temporal decorrelation of turbulent fluctuations, such that it becomes negligibly 
small for T l/?7fc. 

Here we consider Kolmogoroff turbulence for which the energy density spectrum is given by the power law [ssj 

Ek = CKe'/'fc-5/3^ ko<k< fed, (24) 

defined over the range of wavenumbers from fcp, determined by the stirring length scale Ls = 27r/fco on which the 
energy is injected into turbulent motions, to kd, determined by the dissipation length scale Ld = 2ti /kd on which 
the plasma kinetic energy is thermalized. Here Ck is a constant of order unity; for simplicity we set Ck — 1- The 
parameter e is the energy dissipation rate per unit enthalpy, £ ~ Pvac/(''Tw). The corresponding autocorrelation 
function is 1221 



(25) 
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We assume that the stirring and dissipation scales are well separated, i.e., ko <^ k^, which corresponds to the turbulence 
having high Reynolds number. This will be an excellent approximation in any early universe phase transition with the 
stirring scale related to the Hubble length. For simplicity, we adopt Kraichnan's square exponential time dependence 
[3^ to model the temporal decorrelation, 

/fe,r)=cxp(-Jry^2). (26) 

While other forms of /(??fc, t) are also frequently used (see, e.g., [131), neither total power of generated waves nor the 
spectrum are very sensitive to the specific form of the temporal decorrelation [2l| . 

To compute the fourth-order velocity correlation tensors Eq. (jlSp needed in the gravitational wave formula Eq. (I17p , 
we invoke the Millionshchikov quasi-normal hypothesis [2^ : 

{vyylvD = {vy]){vlv1) + (vIvDiv^vt) + (^.^?)(«>^-), (27) 

where vf = Wi(x, t) and = Wi(x + r, t + t). Using Eqs. ([3]), and we obtain 

1 

3^ 

Then Eq. can be evaluated using Eqs. and the convolution theorem to give 

U2 lA h2 1.2 



i?jjij (x',x' + r, r) = Rii{r,T)Rjj{r,T) + -Rij{y,T)R^j{Y,T). (28) 



Hijij{k,u;) = i J dkidujig{ki,uJi)g{k-ki,uj-uji) 



h.2 1.4 h.2 7,2 

27 - _ + — + _ _ + __ 

kl 2klu^ 2^2 m2 ^ 2fc2 



(29) 



where we have defined u — Ik — ki I and 



Ek 



w) = exp ( -— 2 ) . (30) 



Choose the vector k as the axis for spherical coordinates {61, (pi) of the ki integral. The azimuthal angular integral 
over is trivial. The dependence on the direction of ki is clearly only through k • ki, so i/yy (k, cj) = iJ^y (fc, w). 
The other angular integral can be simplified by changing variables from 9i to u, giving 

Hyij{k,uj) = - J dkiduJig{ki,uji) — j J duug{u,Lj - uji) 

+ - dkidu)ig{ki,uJi)[—— + -^-kki] / du -g{u,u; - ui) 
37 V2fci 2fc J J\ki-k\ u 



7T f 1 

+ - J dkidujig{ki,uji)— 



ki+k 

duu^g{u,Lj — ui). (31) 



We need to integrate this expression numerically. The oji integral can be done analytically in terms of the error 
function; the entire expression is reduced to an integral over two dimensionless quantities in Appendix A, Eq. (jA.4[l . 
The result scales with the stirring scale feg, and depends on the Mach number M = {e/koY^^ of the turbulence. Its 
dependence on the dissipation scale kd is through the Reynolds number R — (kd/ko)'^^^; as expected from physical 
considerations, the radiated power is almost completely independent of R. Numerical results are displayed in the next 
Section. 



IV. RELIC GRAVITATIONAL WAVES 



The previous Section and the Appendix has given an analytic expression for the gravitational wave energy spectrum 
resulting from a period of turbulence lasting a time tt, stirred on a scale fco, with Reynolds number R and Mach 
number M. The only significant approximation made is that the turbulence is stationary and acts as a source of 
gravitational waves for a finite time interval; the error made through this idealization is discussed below. In order to 
improve on this approximation, it would be necessary to create a detailed numerical model of the turbulent source, 
including incorporating an actual stirring mechanism, such as colliding bubbles in a phase transition. We have also 
assumed that the expansion of the universe can be ignored during the turbulence; this should be a good approximation 
for any realistic early-universe phase transition. The main effect of expansion would be only to damp the total energy 
in the turbulence by a modest fraction, assuming the turbulence does not last much longer than a Hubble time. 
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A. The Spectrum at the Present Epoch 



To obtain the present spectrum, the gravitational waves generated by the turbulent source must be propagated 
through the expanding universe until today. The wavelengths of the gravitational waves simply scale with the scale 
factor a of the universe, while their total energy density evolves like a~'^ and their amplitude decays like . From 
/OGw(w), Eq. ((2T|) . we can form fJcC'^) = PGwW)l Pc with the critical density pc = 3i/o/87rG. Then, changing to 
linear frequency / = Lojlix^ a characteristic strain amplitude is conventionally defined as 

/ie(/) = 1.263 X 10-18 {hl^G{!)\ (32) 

where is the current Hubble parameter Hq in units of 100 kmsec-^Mpc-^. From the computed hdf) at the epoch 
of the turbulence, given by a scale factor a*, the factor by which the amplitude and the frequency are reduced is 

where is the temperature of the universe with scale factor a*, and is the effective number of relativistic degrees 
of freedom the universe has at this time. To give expressions which are physically transparent, we write the turbulence 
stirring scale and the turbulence duration as fractions of the Hubble length during the turbulence: 

7i/-i = 27r/fco, C^*"' = tt; (34) 

in other words, 7 is the stirring scale's fraction of the Hubble length and C, is the turbulence duration's fraction of the 
Hubble length. For any particular angular frequency of the radiation at the time of the phase transition, we can 
then convert and /ic(w*) to the amplitude hc[f ) and frequency / of the relic gravitational wave background today 
using the useful expressions for a radiation-dominated universe 

4p* 27r^ 4 TT ^ CR 1/2 ^* ,nr\ 

3 45 mpi 

to get 



/ = 1.65 X 10-^ Hz ^ f ^) {^) [ * , (36) 
\koJ\100J VO.Oiy llOOGeV/' ^ ' 



Kif) = 1.28 X 10- (i^) (^^Y' {^f [^y^ [fc3..(/),^.,.,.(.4/),.4/))] . (37) 

The characteristic strain spectrum hdf) is plotted in Fig. 1. The solid lines show three different values for the Mach 
number, M = 0.01, M = 0.1, and M = 1, from lowest to highest amplitude. This dependence on M is in addition to the 
explicit scaling in Eq. (jA.4[) . which is accounted for in the y-axis units. The peak frequency of the spectrum scales 
inversely with the stirring scale and linearly with the characteristic fluid velocity, which is proportional to the Mach 
number. The peak frequency is thus proportional to the inverse of the circulation time on the stirring scale of the 
turbulence. This is the usual result for radiation generation: the characteristic frequency of radiation is determined 
by the characteristic time scale of the source. 

The characteristic parameter values to which the numbers in the plot are scaled (T, = 100 GeV, = 100, 7 — 
C = 0.01) are values consistent with turbulence arising from a strongly first-order phase transition at the electroweak 
scale; see [isj for a detailed discussion of the appropriate parameters. 



B. The Aeroacoustic Limit 



Also plotted in Fig. 1 is an approximation common in aeroacoustics [23|, which replaces Hijij{k = uj,uj) with 
Hijij{k = 0,u}). It is clear that this simplifying approximation is very good for M < 0.1, and overestimates the 
maximum amplitude of hc{f) by around 30% for M = 1. In this limit, the argument of the Bessel function in Eq. (PD|) 
becomes small. Substituting Eq. (PU]) into Eq. (PT|) gives 

PGwico) = 4c^3gwVT / dTdU^e''^^joiujORmi^,T). (38) 
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FIG. 1: The spectrum of gravitational radiation from turbulence. The three solid lines are for different Mach numbers, with 
M = 0.01, M = 0.1, and M = 1 from lowest to highest amplitude. Note that these three cases have also been scaled by a factor of 
M~^^'^ for display, since this is how the low-frequency tail scales with M. The dotted lines, which are virtually indistinguishable 
from the solid lines except for the M = 1 case, show the fc = approximation to the gravitational wave source. 



Thus if uj£, is small compared to unity, the aeroacoustic limit fc ^ is guaranteed to be valid. 

In the case of aeroacoustics, this approximation works because the fluid velocity is always assumed to be small 
compared to the velocity of the radiated acoustic waves (low Mach number). In the cosmological regime, the interesting 
case is for plasma with a relativistic amount of kinetic energy (otherwise there is not substantial gravitational radiation 
produced). This will occur only when the plasma is at a high enough temperature that it is fully relativistic: otherwise, 
the amount of energy injected into plasma motions would have to be a substantial fraction of the particle mass scale 
rather than the cosmological temperature scale, and this is unlikely on general grounds. A relativistic plasma has 
sound speed 1/V5, and the Mach number of the turbulent plasma can never be much larger than 1; it will also not 
be too much smaller than 1. In this case, the fluid velocities will be roughly the sound speed, but this is close to 
the propagation speed of the emitted radiation. Therefore, we do not automatically have ^cj <C 1 in Eq. (|38p . and 
the validity of the aeroacoustic approximation must be ascertained by explicit calculation. As we see in Fig. 1, the 



8 



approximation still gives the right order of magnitude for the spectrum amplitude even for Mach number M = 1, 
corresponding to a fluid velocity equal to the sound speed. 



C. Asymptotic Limits 

The validity of the aeroacoustic approximation /c — > simplifles finding asymptotic forms for the spectrum. Consider 
various frequency regimes of Eq. ([3T|) with fc = 0; this is facilitated by Eq. (jA.sp . Note that in this limit, the dependence 
on the Mach number M simply scales with the frequency. We assume R ^ 1, or else fully developed turbulence cannot 
exist; this is an excellent approximation for early-universe plasma stirred on scales near the Hubble length. In the 
low-frequency regime, simply take the limit a; ^ to get 



7M3 

60fc4(7)3/2 



(0,0.)- ^^,^4,_,,/2 ' ^ ^ 0. (39) 



Physically, these frequencies are lower than the lowest characteristic frequency in the problem, corresponding to the 
eddy turnover time on the stirring scale. This result of a constant Hijij is universal and does not depend on either 
the spectrum or temporal characteristics of the turbulence (see, e.g., Refs. [13, [3 [3 [3 ^'^'^ ^^^^ IS)- It translates 
to hc{f) scaling as /^/^ at low /. 

At high frequencies w ^ fcoMR^/^, the integral in Eq. (jA.5[) is dominated by the contribution from its lower limit. 
After using the asymptotic form erfc(a;) ~ x~^it~^/'^ exp(— a;^), x — > oo, an integration by parts gives the leading-order 
asymptotic behavior as 



■7\/\i h^M^ 

''"'^^''^^ ^ 32^2fc„4RV4 1,2 exp(-2c.V(fcgM2R)) , u » fcoMR^^. (40) 

This exponential suppression is evident in Fig. 1; the dependence on R is negligible as expected from physical 
considerations. The functional form of the high-frequency suppression is determined by the specific form of the time 
autocorrelation function of the turbulence, Eq. (|26p . but for any autocorrelation the amplitude of the emitted waves 
should be very small in this regime. Physically, this limit corresponds to radiation frequencies which are larger than 
any frequencies in the turbulent motions; consequently, no scale of turbulent fluctuations generates these radiation 
frequencies directly, and the resulting small radiation amplitude is due to the sum of small contributions from many 
lower-frequency source modes. Since the integral is dominated by the lower integration limit, the highest-frequency 
source fluctuations (which contain very little of the total turbulent energy) contribute most to this high-frequency 
radiation tail. 

In the intermediate frequency regime, for frequencies fcoM < a; < fcoMR^/^, the integral in Eq. llA.SP is dominated 
by the contribution around fepM^/w^ due to the exponential factor in the integrand, with a width of the same order. 
Physically, this implies that radiation emission at some frequency in this range is dominated by the turbulent vortices 
of the same frequency. Consequently, we have the rough estimate 



This yields /ic(/) oc /~^^^'', compared to /ic(/) oc f^^lJ"^ f^^"^/^ in Ref. jL4.], where /stir is the turbulence circulation 
frequency at the stirring scale. The slight discrepancy from the spectrum shape in Ref. [13] comes about because we 
have treated the time correlations of the turbulence in a more realistic way. Here we distinguish two time scales, the 
decorrelation time which describes how the fluid velocities in a given size eddy are correlated with each other after a 
given time interval, and the largest eddy turnover time. In practice, the dropoff with frequency in this regime is strong 
enough that the high-frequency behavior in Eq. (1401) only holds when the spectrum is many orders of magnitude below 
the peak amplitude. 

The intermediate frequency regime scaling with frequency depends on the specific model of the turbulence power 
spectrum. The Kolmogoroff model is not the only possibility, especially in the presence of magnetic fields. Any 
model of turbulence which includes the local transfer of energy in the wavenumber space will satisfy El/rjk oc fc~'^ 
[soj . In the fc = limit, it is straightforward to derive that in general, Hijij(0,q) scales as l/g^^", where n is the 
scaling exponent of the turbulence autocorrelation function. For Kolmogoroff turbulence, n = 2/3 (Eq. ([25]) ). But 
for Iroshnikov-Kraichnan turbulence (for example), n = 1, and consequently the frequency dependence is somewhat 
softer, Hijij{0, q) oc In practical terms, this modified turbulence spectrum produces radiation with very similar 

detectability properties to that from the Kolmogoroff turbulence spectrum considered here. 
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V. DISCUSSION 

We have calculated for the first time the spectrum of relic gravitational radiation resulting from a period of stationary 
turbulence in the early universe, in terms of the turbulence duration, stirring scale, Reynolds and Mach numbers, 
and the temperature of the universe when the turbulence occurs. This is probably the best that can be done without 
a detailed simulation of actual turbulent motions. The most likely source of energy injection leading to turbulence 
is an early-universe phase transition; the connection between the phenomenological parameters describing a phase 
transition and the parameters describing the turbulence are given explicitly in Ref. [14i] . 

The calculation we present here is conceptually simple. The only assumptions made are that the turbulence lasts 
for a finite duration which is at least a turnover time on the stirring scale, and that during this time the turbulence 
can be characterized as stationary. The spatial power spectrum is taken to be the Kolmogoroff form, Eq. (|24|) . with 
temporal correlations of the Kraichnan form, Eq. ([26|) . These scalings are appropriate for non-relativistic turbulence 
with large Reynolds number. While the cosmological case will have large Reynolds numbers, the turbulence will be 
relativistic in the most interesting cases for gravitational radiation generation. As argued in Ref. [l^ . a non-relativistic 
approximation to relativistic turbulence likely underestimates the resulting radiation: relativistic turbulence contains 
more kinetic energy for a given fiuid velocity. We expect the same general results to hold, except the expression for the 
Mach number M'^ = e/fcp will clearly be modified, giving smaller Mach numbers than this non-relativistic expression. 

Cosmological turbulence will never be precisely stationary, since the universe is expanding. Turbulence from a 
phase transition will also not be stationary because the duration of the source is comparable to the eddy turnover 
time on the stirring scale so the turbulence will decay with time. Even so, as long as the eddies on a given length 
scale can be treated as uncorrelated sources of turbulence, Ref. [l^ argues that the resulting radiation spectrum will 
be close to that from a stationary source, simply due to the inevitable cascade of energy from the stirring scale down 
to the diffusion scale. This point can be made somewhat more formally, using an argument similar to that given 
by Proudman [U, In the case of stationary turbulence, the time derivatives in Eq. ([7]) lead to factors of I/tq 
when computing the radiation spectrum. If the turbulence is decaying, then additional terms proportional to time 
derivatives of the correlation functions also will appear. But the characteristic time scale of the turbulence decay 
is at least several times greater than the turnover time on the stirring scale, and consequently, these additional terms 
which are proportional to l/r^ can be neglected compared to the stationary term. 

It has been claimed that turbulence in the early universe is a source of such short duration that it cannot be treated 
as stationary at all, and that the resulting radiation spectrum should be imprinted with the characteristic wave 
number of the turbulent source instead of its characteristic frequency [111 . We fully agree that a gravitational wave 
source lasting for a sufficiently short duration will not be well described by a short piece of a stationary source; clearly 
in the limit that the turbulence duration r goes to zero, the resulting radiation spectrum should peak at a frequency 
corresponding to the characteristic wavenumber of the source (the inverse stirring scale) , and our formalism will not be 
valid in that limit. However, it is straightforward to see that a turbulent source in the early universe will actually last 
long enough so that our calculations are valid. Ref. [sl] uses a simple toy model for a stochastic cosmological gravity 
wave source to argue that the condition ttlOs w 1 represents the dividing line between sources that imprint their 
characteristic frequency on the radiation and could be described using our formalism, and sources that imprint their 
wavenumber; here we write the duration of the source as tt and the characteristic source angular frequency as ujs- If 
we write the turbulence turnover time on the stirring scale as tstir, the associated angular frequency is lOs = 27r/tstir; 
also write the duration of the turbulent source as tt = Ntstir, so that it lasts N times the stirring-scale turnover 
time. We can thus write lOsTt — 2TrN, and even for the unrealistically short duration = 1, uJsTt is significantly 
larger than unity. For realistic turbulence, we expect the dissipation time to be multiple turnover times. Inspecting 
the exact gravitational wave solution for the toy model source in Ref. |38 | confirms that, even for N as small as 1, 
the radiation spectrum will be peaked at tUs- 

We also assume that turbulence is non-magnetic and non-helical. Either of these complications can modify the 
power law in Eq. ([M)) or the form of the time correlation Eq. (pS)) j2^. The main effect of any modification is 
to change the rate at which the radiation spectrum falls off at high frequencies, but since this dependence is quite 
steep, even substantial changes to the asymptotic behavior of the spectrum lead to little qualitative difference in the 
spectrum. As mentioned above, the low-frequency behavior is independent of any details of the turbulence, and the 
peak frequency is determined by the eddy turnover time on the stirring scale where the energy density peaks, which 
will also be independent of any details of the turbulent cascade. 

The proposed Laser Interferometer Space Antenna (LISA) satellite mission has a 5ct strain sensitivity to stochastic 
backgrounds of below he = 10~^^ between frequencies 10"'^ and 10~^ Hz, and decreasing to around he ~ 10~^° at 
10"'* Hz, for one year of integration (see, e.g., [40|). Comparing with Fig. 1, turbulence with a Mach number M = 1 
would be a factor of 1000 larger than the LISA detection threshold at the peak frequency around 10^'' Hz. For a 
Mach number M = 0.1, the peak amplitude decreases by a factor of 100 due to the M~'^/^ scaling and the different 
signal spectrum. However, the peak frequency also shifts to 10~^ Hz, at which point LISA's sensitivity has declined 
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greatly; the steep high-frequency tail of the gravitational wave spectrum makes detection with LISA marginal in 
this case. Detectors consisting of two or more correlated LISA detectors or enhanced versions of LISA optimized 
for detecting stochastic backgrounds have been discussed [4l[, such as the envisioned GREAT mission [42]; future 
space-based interferometers could be configured to give strain sensitivities comparable to LISA, but with a frequency 
window between 10~* and 10^® Hz. Such an experiment would easily detect turbulence at the electroweak scale with 
a Mach number M = 0.1, and would even flirt with a detection at M = 0.01. Turbulence generated at somewhat 
higher energy scales shifts to higher frequencies and easier detection with LISA. 

As is widely appreciated, detecting cosmological backgrounds of gravitational radiation is not only an issue of 
detector sensitivity, but also of foreground discrimination. The galactic population of short-period binaries of compact 
objects, mostly white dwarfs, is known to produce a confusion-limited stochastic background at frequencies below 10"'^ 
Hz [4^. At low frequencies, separating this galactic source from a cosmological source is essential, likely by exploiting 
the non- uniform directional distribution of an galactic source 0, . A uniform stochastic source arising from the 

confusion limit of numerous extragalactic binaries provides a further complication (47j , which can only be distinguished 
from a primordial background via differing spectra. We also note that the source of the turbulence itself may produce 
a gravitational wave spectrum, and that the characteristic peak frequency may scale differently from the turbulent 
spectrum; see, e.g., the spectra for first-order phase transitions in Ref. [13|. A distinctive two-peaked shape to the 
gravitational wave spectrum in certain regions of parameter space will also aid in its detection. 

We have no guarantees of violent events in the early universe. However, turbulence is a completely generic result of 
energy injection on a characteristic length scale, and we have shown in this paper that the resulting relic gravitational 
waves are within the realm of detectability, even for turbulence with Mach numbers as low as 0.01, corresponding to 
an energy input into the early universe of 10"'* of the total energy density. Many scenarios for the electroweak phase 
transition [i^ and other physics [2^ will result in releases of energy that are interestingly large. The remarkable 
possibility of probing high-energy physics via the detection of vanishingly small spacetime distortions left from when 
the universe was a trillionth of a second old impels us to look. 



APPENDIX: NUMERICAL EVALUATION OF Hijij{k,uj) 

We need to evaluate Eq. pT|) explicitly, with g{k,uj) given by Eq. ([30| . The integral over uji can be evaluated 
analytically using the identity 



f°° 1 
J dy expi-Ay"^) exp{~B{x - yf) = - 



2\A + B 



1/2 



exp 



A + B 



erfc — 



Bx 



VA + B 



(A.l) 



This expression is simple to derive by writing the integrand as a single exponential and completing the square in 
the argument of the exponential, followed by a linear change of variables to give the error function. Then Eq. (|3ip 
becomes 



247r(27r)3/2fc 



X exp 



I 



-1/2 



27- 



fc2 



2^2 



2£-2/3a;2 



erfc 



2l/2£-l/3, 



4/3 



-4/3 



1/2 



2fc2 
(A.2) 



The lower limit on the u integral is niax[|A:i — fc|, /cq] and the upper limit is niin[fci + fc, /c^], provided the lower limit is 
less than the upper limit; otherwise the integral over u is zero. These conditions on the limits arise due to the limited 
range of fc over which the function E-k has support. Note that Eq. (|A.2p is regular as fc — > 0, with the limit 



y (0, W) 



Is 



dkik^ ^ exp 



,2/3fc4/3^ 



erfc 



^l/3fc2/3^ 



(A.3) 



Now rescale all dimensionful quantities by powers of ka to make them diniensionless; we abbreviate e/fco = M'^, 
where M is the Mach number of the turbulence, fcd/fco = R'^^'*, where R is the Reynolds number of the turbulence. 



p = k/ko, and q = Lu/ko- The change of variables x = (fci/fcc 



\-4/3 



(u/ko) 4/3 simplifies the remaining integrals. 



giving 



Hij.j (p, q) 



SM^fcn 



256(27r)3/2p 



dxx^^^ / dyy^/^ix + y)-^/"^" 



exp 



2xy g2 
x -f J/ M2 



erfc — 



21/2, 
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54 - 2p2a;3/2 _ 2p2y3/2 ^ p4^3/2^3/2 + ±_ + 



^3/2 
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the lower limit of the y integral is max[(a; 



-3/4 



-4/3 



R ^] and the upper limit is min[ 



-3/4 



-4/3 



(A.4) 



1], provided 



the lower limit is less than the upper limit. In the limit p — > 0, both of these limits are x, so the integral has a leading 
order behavior proportional to p and thus Hijij (p, q) is regular, with the limit 



Hzjij (0, q) 



167r3/2 



dx x^^^'^ exp (— g^x) erfc (^—qx^^^^ 



(A.5) 



where we have abbreviated q = q/M since in this limit the integral depends only on q and not on either g or M 
separately, aside from the constant prefactor. Note that R 3> 1 for a medium which supports turbulence; we expect 
R > 2000 during the cosmological epochs of relevance. The integrals converge as i? ^ 00, and the lower limit of the 
x-integrals in Eqs. (jA.4p and (jA.5[) can be replaced by zero. In Eq. (|A.4[) . the terms with factors of x^^^^ and y^^^^ 
in the integrand converge somewhat slowly but have small prefactors compared to the first term, giving a negligible 
dependence of the integral on the diffusion scale. Numerically, it is convenient to take R as some large but finite 
value; then the integrand in Eq. (jA.4[) is smooth and regular over the full range of integration, and can now be easily 
performed for any values of p and q. 
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